{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flexural Subsidence\n",
    "\n",
    "This example explores how to use a BMI implementation using sedflux's subsidence model as an example.\n",
    "\n",
    "### Links\n",
    "* [sedflux source code](https://github.com/mcflugen/sedflux): Look at the files that have *deltas* in their name.\n",
    "* [sedflux description on CSDMS](https://csdms.colorado.edu/wiki/Model_help:Sedflux): Detailed information on the CEM model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interacting with the Subside BMI using Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some magic that allows us to view images within the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the `Subside` class, and instantiate it. In Python, a model with a BMI will have no arguments for its constructor. Note that although the class has been instantiated, it's not yet ready to be run. We'll get to that later!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymt.models\n",
    "subside = pymt.models.Subside()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though we can't run our subsidence model yet, we can still get some information about it. *Just don't try to run it.* Some things we can do with our model are get the names of the input variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.output_var_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.input_var_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get information about specific variables. Here we'll look at some info about lithospheric deflections. This is the main input of the Subside model. Notice that BMI components always use [CSDMS standard names](http://csdms.colorado.edu/wiki/CSDMS_Standard_Names). With that name we can get information about that variable and the grid that it is on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK. We're finally ready to run the model. Well not quite. First we initialize the model with the BMI **initialize** method. Normally we would pass it a string that represents the name of an input file. For this example we'll pass **None**, which tells Cem to use some defaults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "config_file, config_folder = subside.setup()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.initialize(config_file, dir=config_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.var[\"earth_material_load__pressure\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.grid[0].node_shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running the model, let's set an input parameter - the overlying load."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load = np.zeros((500, 500))\n",
    "load[250, 250] = 1e3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main output variable for this model is *deflection*. In this case, the CSDMS Standard Name is:\n",
    "\n",
    "    \"lithosphere__increment_of_elevation\"\n",
    "\n",
    "First we find out which of Subside's grids contains deflection. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.var['lithosphere__increment_of_elevation']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the *grid_id*, we can now get information about the grid. For instance, the number of dimension and the type of grid (structured, unstructured, etc.). This grid happens to be *uniform rectilinear*. If you were to look at the \"grid\" types for wave height and period, you would see that they aren't on grids at all but instead are scalars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.grid[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because this grid is uniform rectilinear, it is described by a set of BMI methods that are only available for grids of this type. These methods include:\n",
    "* get_grid_shape\n",
    "* get_grid_spacing\n",
    "* get_grid_origin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allocate memory for the water depth grid and get the current values from `cem`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.set_value(\"earth_material_load__pressure\", load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subside.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dz = subside.get_value('lithosphere__increment_of_elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(dz.reshape((500, 500)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load[125, 125] = 2e3\n",
    "subside.set_value(\"earth_material_load__pressure\", load)\n",
    "subside.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dz = subside.get_value('lithosphere__increment_of_elevation')\n",
    "plt.imshow(dz.reshape((500, 500)))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
